Practice Dec 12, 2016
Problem 1
Remember the tomato data set generated by Pepe; we first looked at this when we were working on Chapter 9. 35 accessions for seven species were grown in sun and shade.
Assess whether there is evidence for total length (“totleng”) response to shade and whether this response varies by species. Consider whether including accession (“acs”) or shelf (“shelf”) using adaptive priors improves the model fit.
get the data and format
tomato <- read.csv("Assignment_Chapter_09/TomatoR2CSHL.csv")
names(tomato)## [1] "shelf" "flat" "col" "row" "acs" "trt"
## [7] "days" "date" "hyp" "int1" "int2" "int3"
## [13] "int4" "intleng" "totleng" "petleng" "leafleng" "leafwid"
## [19] "leafnum" "ndvi" "lat" "lon" "alt" "species"
## [25] "who"tomato$species_id <- coerce_index(tomato$species)
tomato$shade <- ifelse(tomato$trt=="L",1,0)
tomato$shelf_id <- coerce_index(tomato$shelf)
tomato$species_id <- coerce_index(tomato$species)
tomato$acs_id <- coerce_index(tomato$acs)
tomato.small <- tomato[,c("shelf_id","acs_id","shade","totleng","species_id")]
summary(tomato.small)## shelf_id acs_id shade totleng
## Min. :1.000 Min. : 1.0 Min. :0.0000 Min. : 13.59
## 1st Qu.:2.000 1st Qu.: 9.0 1st Qu.:0.0000 1st Qu.: 39.25
## Median :3.000 Median :18.0 Median :1.0000 Median : 50.98
## Mean :3.512 Mean :18.1 Mean :0.5089 Mean : 53.70
## 3rd Qu.:5.000 3rd Qu.:27.0 3rd Qu.:1.0000 3rd Qu.: 64.76
## Max. :6.000 Max. :36.0 Max. :1.0000 Max. :129.43
## species_id
## Min. :1.000
## 1st Qu.:2.000
## Median :3.000
## Mean :2.927
## 3rd Qu.:4.000
## Max. :5.000shade effect
mean(tomato.small$totleng[tomato.small$shade==0])## [1] 42.0596m.shade <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
beta_shade*shade,
alpha ~ dnorm(42,10),
beta_shade ~ dnorm(0,20),
sigma ~ dcauchy(0,1)
),
data = tomato.small,
chains=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.364957 seconds (Warm-up)
## 0.304377 seconds (Sampling)
## 0.669334 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.4212 seconds (Warm-up)
## 0.293224 seconds (Sampling)
## 0.714424 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.300826 seconds (Warm-up)
## 0.244408 seconds (Sampling)
## 0.545234 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.274527 seconds (Warm-up)
## 0.257516 seconds (Sampling)
## 0.532043 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000122 seconds (Sampling)
## 0.000126 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]plot(m.shade)
pairs(m.shade)
precis(m.shade)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 42.08 0.69 40.95 43.15 1803 1
## beta_shade 22.81 0.95 21.30 24.33 1862 1
## sigma 16.05 0.34 15.51 16.59 2627 1shade and species
m.shade.species <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
beta_sp[species_id] +
beta_shade*shade,
alpha ~ dnorm(42,10),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
sigma ~ dcauchy(0,1)
),
data = tomato.small,
chains=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 2.65907 seconds (Warm-up)
## 1.81048 seconds (Sampling)
## 4.46955 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 2.69747 seconds (Warm-up)
## 1.57889 seconds (Sampling)
## 4.27636 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 2.71582 seconds (Warm-up)
## 1.84126 seconds (Sampling)
## 4.55708 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 2.66713 seconds (Warm-up)
## 1.65344 seconds (Sampling)
## 4.32057 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 0.000182 seconds (Sampling)
## 0.000187 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]plot(m.shade.species)
pairs(m.shade.species)
precis(m.shade.species,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.13 4.21 34.38 47.85 641 1.01
## beta_sp[1] -1.29 4.25 -8.31 5.34 444 1.00
## beta_sp[2] 0.12 4.27 -7.08 6.52 659 1.00
## beta_sp[3] 2.42 4.25 -4.12 9.43 663 1.00
## beta_sp[4] -9.41 4.29 -16.71 -2.86 675 1.00
## beta_sp[5] 8.35 4.24 1.44 15.05 667 1.00
## beta_shade 22.99 0.93 21.56 24.48 1757 1.00
## sigma 15.17 0.33 14.66 15.69 1947 1.00compare(m.shade,m.shade.species) ## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species 8348.1 6.9 0.0 1 53.11 NA
## m.shade 8458.4 3.2 110.3 0 52.47 21.97levels(tomato$species)## [1] "S. chilense" "S. chmielewskii" "S. habrochaites" "S. pennellii"
## [5] "S. peruvianum"Good evidence for species being important
consider species X shade interaction
m.shade.species.int <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
beta_sp[species_id] +
beta_shade*shade +
beta_sp_shade[species_id]*shade,
alpha ~ dnorm(53,20),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
beta_sp_shade[species_id] ~ dnorm(0,10),
sigma ~ dcauchy(0,1)
),
data = tomato.small,
chains=4,
cores=2
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 8e-06 seconds (Warm-up)
## 0.000484 seconds (Sampling)
## 0.000492 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]plot(m.shade.species.int)
par(mfrow=c(1,1))
precis(m.shade.species.int,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.58 4.42 34.51 48.32 990 1
## beta_sp[1] -4.24 4.56 -11.35 2.86 1049 1
## beta_sp[2] 2.67 4.54 -3.99 10.04 1044 1
## beta_sp[3] 3.15 4.56 -4.02 10.11 1040 1
## beta_sp[4] -11.83 4.61 -19.03 -4.53 1058 1
## beta_sp[5] 7.03 4.55 -0.15 14.12 1081 1
## beta_shade 22.19 4.27 15.63 29.17 1075 1
## beta_sp_shade[1] 5.84 4.63 -1.23 13.56 1219 1
## beta_sp_shade[2] -5.10 4.47 -12.07 2.09 1224 1
## beta_sp_shade[3] -1.71 4.53 -9.20 5.49 1180 1
## beta_sp_shade[4] 4.43 4.68 -2.83 12.31 1344 1
## beta_sp_shade[5] 2.51 4.57 -4.65 9.80 1194 1
## sigma 15.03 0.33 14.51 15.53 2566 1levels(tomato$species)## [1] "S. chilense" "S. chmielewskii" "S. habrochaites" "S. pennellii"
## [5] "S. peruvianum"compare(m.shade,m.shade.species,m.shade.species.int)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.int 8334.6 10.6 0.0 1 53.30 NA
## m.shade.species 8348.1 6.9 13.5 0 53.11 8.16
## m.shade 8458.4 3.2 123.8 0 52.47 22.14plot(coeftab(m.shade,m.shade.species,m.shade.species.int))
Reasonable evidence for the interaction, although the 95% posterior intervals overlap among all of the species_shade interactions
Add adaptive prior for accessions
m.shade.species.int.acs <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
alpha_acs[acs_id] +
beta_sp[species_id] +
beta_shade*shade +
beta_sp_shade[species_id]*shade,
alpha ~ dnorm(53,20),
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
beta_sp_shade[species_id] ~ dnorm(0,10),
c(sigma,sigma_acs) ~ dcauchy(0,1)
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=2
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.000359 seconds (Sampling)
## 0.000362 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]plot(m.shade.species.int.acs,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
par(mfrow=c(1,1))
precis(m.shade.species.int.acs,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.50 4.44 34.56 48.78 2782 1
## alpha_acs[1] -5.40 3.68 -11.19 0.52 4841 1
## alpha_acs[2] -10.63 3.53 -16.48 -5.21 4319 1
## alpha_acs[3] -3.54 3.42 -9.01 1.83 4277 1
## alpha_acs[4] 3.15 3.38 -1.95 8.75 4305 1
## alpha_acs[5] 13.56 3.41 8.13 18.93 4221 1
## alpha_acs[6] 2.48 3.44 -2.97 7.87 4133 1
## alpha_acs[7] 9.86 3.46 4.30 15.31 4518 1
## alpha_acs[8] 6.35 3.33 1.24 11.84 3916 1
## alpha_acs[9] 3.44 3.37 -1.59 9.08 3752 1
## alpha_acs[10] 3.86 3.90 -2.12 10.29 5649 1
## alpha_acs[11] 2.31 3.32 -2.90 7.74 3848 1
## alpha_acs[12] -9.78 3.50 -15.44 -4.35 4320 1
## alpha_acs[13] -4.01 3.66 -9.69 1.93 4923 1
## alpha_acs[14] -0.69 3.53 -6.04 5.22 4108 1
## alpha_acs[15] -3.50 4.99 -11.56 4.29 7840 1
## alpha_acs[16] 1.36 3.29 -3.76 6.70 4028 1
## alpha_acs[17] 4.78 3.39 -0.61 10.21 4564 1
## alpha_acs[18] -12.07 3.53 -17.48 -6.29 4835 1
## alpha_acs[19] 15.04 3.46 9.71 20.77 4066 1
## alpha_acs[20] -6.18 3.31 -11.25 -0.70 3712 1
## alpha_acs[21] 0.96 3.36 -4.43 6.34 3990 1
## alpha_acs[22] 6.18 4.52 -1.01 13.33 6243 1
## alpha_acs[23] -2.24 3.68 -8.12 3.56 4884 1
## alpha_acs[24] -4.28 3.38 -9.49 1.29 4363 1
## alpha_acs[25] -5.10 3.31 -10.27 0.32 3994 1
## alpha_acs[26] -0.31 3.55 -5.93 5.42 4978 1
## alpha_acs[27] 1.83 3.33 -3.51 7.01 4483 1
## alpha_acs[28] -4.60 3.47 -10.44 0.68 4768 1
## alpha_acs[29] 1.14 3.55 -4.56 6.74 4772 1
## alpha_acs[30] 6.78 3.42 1.31 12.18 4699 1
## alpha_acs[31] 10.34 3.49 4.88 15.93 3991 1
## alpha_acs[32] 0.51 3.59 -5.26 6.08 4599 1
## alpha_acs[33] -3.30 3.54 -8.88 2.39 4628 1
## alpha_acs[34] -8.61 3.40 -14.11 -3.28 4360 1
## alpha_acs[35] -5.57 3.53 -11.14 -0.03 4579 1
## alpha_acs[36] -7.25 4.36 -14.21 -0.35 5457 1
## beta_sp[1] -4.01 5.00 -12.29 3.59 2889 1
## beta_sp[2] 2.54 4.98 -5.24 10.57 2864 1
## beta_sp[3] 2.65 4.93 -5.24 10.52 2974 1
## beta_sp[4] -10.69 5.09 -18.66 -2.42 3197 1
## beta_sp[5] 6.23 5.04 -2.03 14.05 2759 1
## beta_shade 22.05 4.43 14.73 28.86 2419 1
## beta_sp_shade[1] 5.72 4.66 -2.02 12.86 2695 1
## beta_sp_shade[2] -5.12 4.64 -11.88 2.93 2648 1
## beta_sp_shade[3] -1.89 4.65 -8.95 5.86 2676 1
## beta_sp_shade[4] 5.11 4.80 -2.46 12.78 2827 1
## beta_sp_shade[5] 2.39 4.65 -4.59 10.31 2744 1
## sigma 13.34 0.31 12.85 13.82 8622 1
## sigma_acs 7.47 1.08 5.70 9.05 5679 1compare(m.shade.species.int,m.shade.species.int.acs)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.int.acs 8124.3 38.7 0.0 1 55.22 NA
## m.shade.species.int 8334.6 10.6 210.3 0 53.30 28.15coeftab(m.shade.species.int,m.shade.species.int.acs)## m.shade.species.int m.shade.species.int.acs
## alpha 41.58 41.50
## beta_sp[1] -4.24 -4.01
## beta_sp[2] 2.67 2.54
## beta_sp[3] 3.15 2.65
## beta_sp[4] -11.83 -10.69
## beta_sp[5] 7.03 6.23
## beta_shade 22.19 22.05
## beta_sp_shade[1] 5.84 5.72
## beta_sp_shade[2] -5.10 -5.12
## beta_sp_shade[3] -1.71 -1.89
## beta_sp_shade[4] 4.43 5.11
## beta_sp_shade[5] 2.51 2.39
## sigma 15.03 13.34
## alpha_acs[1] NA -5.4
## alpha_acs[2] NA -10.63
## alpha_acs[3] NA -3.54
## alpha_acs[4] NA 3.15
## alpha_acs[5] NA 13.56
## alpha_acs[6] NA 2.48
## alpha_acs[7] NA 9.86
## alpha_acs[8] NA 6.35
## alpha_acs[9] NA 3.44
## alpha_acs[10] NA 3.86
## alpha_acs[11] NA 2.31
## alpha_acs[12] NA -9.78
## alpha_acs[13] NA -4.01
## alpha_acs[14] NA -0.69
## alpha_acs[15] NA -3.5
## alpha_acs[16] NA 1.36
## alpha_acs[17] NA 4.78
## alpha_acs[18] NA -12.07
## alpha_acs[19] NA 15.04
## alpha_acs[20] NA -6.18
## alpha_acs[21] NA 0.96
## alpha_acs[22] NA 6.18
## alpha_acs[23] NA -2.24
## alpha_acs[24] NA -4.28
## alpha_acs[25] NA -5.1
## alpha_acs[26] NA -0.31
## alpha_acs[27] NA 1.83
## alpha_acs[28] NA -4.6
## alpha_acs[29] NA 1.14
## alpha_acs[30] NA 6.78
## alpha_acs[31] NA 10.34
## alpha_acs[32] NA 0.51
## alpha_acs[33] NA -3.3
## alpha_acs[34] NA -8.61
## alpha_acs[35] NA -5.57
## alpha_acs[36] NA -7.25
## sigma_acs NA 7.47
## nobs 1008 1008plot(coeftab(m.shade.species.int,m.shade.species.int.acs),pars=1:13)
Adding acs does not change most of the shared coefficients by much, although it drops sigma.
Add shelf random effect
table(tomato.small$shade, tomato.small$shelf_id)##
## 1 2 3 4 5 6
## 0 0 0 0 174 125 196
## 1 161 174 178 0 0 0table(tomato.small$acs_id,tomato.small$shelf_id)##
## 1 2 3 4 5 6
## 1 0 5 5 0 4 6
## 2 2 6 6 3 2 6
## 3 6 6 6 6 4 6
## 4 6 6 6 6 3 6
## 5 6 6 6 5 5 6
## 6 5 5 5 5 5 6
## 7 5 6 6 6 2 6
## 8 7 4 6 9 5 6
## 9 6 5 6 6 5 6
## 10 0 6 2 0 3 5
## 11 6 6 6 9 4 5
## 12 2 4 6 3 5 6
## 13 3 6 5 1 3 6
## 14 2 2 6 3 4 6
## 15 0 2 0 0 0 3
## 16 8 6 6 9 5 6
## 17 6 5 5 7 4 6
## 18 4 6 3 5 4 6
## 19 5 6 5 5 3 6
## 20 10 2 6 11 3 5
## 21 8 4 4 9 5 6
## 22 0 5 0 0 0 3
## 23 3 4 5 0 3 6
## 24 6 6 6 6 5 6
## 25 8 6 6 9 4 6
## 26 3 3 6 4 2 4
## 27 7 6 5 9 4 6
## 28 4 6 5 5 4 5
## 29 4 4 5 5 3 5
## 30 5 6 7 5 3 6
## 31 6 3 6 4 4 5
## 32 3 4 4 6 3 5
## 33 5 5 6 5 1 6
## 34 5 6 5 5 4 6
## 35 3 5 3 3 4 6
## 36 2 1 3 0 3 1most accessions are on most shelves
m.shade.species.int.acs.shelf <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_shade*shade +
beta_sp_shade[species_id]*shade,
alpha ~ dnorm(53,20),
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
beta_sp_shade[species_id] ~ dnorm(0,10),
c(sigma,sigma_acs) ~ dcauchy(0,1),
sigma_shelf ~ dexp(1) # with only 3 shelves per trt might be hard to estimate w. cauchy
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=2
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 7e-06 seconds (Warm-up)
## 0.000508 seconds (Sampling)
## 0.000515 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image(file="tomato1212.Rdata")plot(m.shade.species.int.acs.shelf,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
par(mfrow=c(1,1))
precis(m.shade.species.int.acs.shelf,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.28 4.71 33.63 48.64 3876 1
## alpha_acs[1] -6.46 3.66 -12.32 -0.60 6404 1
## alpha_acs[2] -11.35 3.56 -16.99 -5.65 5773 1
## alpha_acs[3] -3.32 3.37 -8.74 1.94 5714 1
## alpha_acs[4] 3.43 3.41 -2.09 8.70 6028 1
## alpha_acs[5] 13.88 3.42 8.44 19.30 5984 1
## alpha_acs[6] 2.73 3.41 -2.48 8.33 6098 1
## alpha_acs[7] 9.93 3.44 4.42 15.36 5701 1
## alpha_acs[8] 6.89 3.30 1.68 12.23 4995 1
## alpha_acs[9] 3.63 3.34 -1.91 8.76 5073 1
## alpha_acs[10] 2.97 3.91 -2.85 9.52 7398 1
## alpha_acs[11] 2.70 3.34 -2.67 7.97 5772 1
## alpha_acs[12] -10.18 3.47 -15.85 -4.70 5454 1
## alpha_acs[13] -4.09 3.66 -9.86 1.74 6426 1
## alpha_acs[14] -0.79 3.55 -6.38 4.90 5565 1
## alpha_acs[15] -5.53 4.92 -12.99 2.72 7819 1
## alpha_acs[16] 1.86 3.33 -3.30 7.17 5622 1
## alpha_acs[17] 5.00 3.43 -0.51 10.43 5559 1
## alpha_acs[18] -12.38 3.55 -17.90 -6.55 5381 1
## alpha_acs[19] 15.11 3.48 9.61 20.59 5765 1
## alpha_acs[20] -5.13 3.32 -10.55 0.06 5262 1
## alpha_acs[21] 1.60 3.30 -3.79 6.74 5078 1
## alpha_acs[22] 4.56 4.52 -2.25 12.24 9067 1
## alpha_acs[23] -2.21 3.74 -8.03 3.90 6604 1
## alpha_acs[24] -4.02 3.35 -9.33 1.34 5907 1
## alpha_acs[25] -4.56 3.30 -9.54 0.92 5635 1
## alpha_acs[26] -0.19 3.63 -5.65 6.00 5843 1
## alpha_acs[27] 2.19 3.38 -2.99 7.75 5346 1
## alpha_acs[28] -4.73 3.48 -10.31 0.86 5552 1
## alpha_acs[29] 1.26 3.55 -4.22 7.07 5770 1
## alpha_acs[30] 6.76 3.44 1.18 12.09 5769 1
## alpha_acs[31] 10.76 3.48 5.45 16.47 5634 1
## alpha_acs[32] 1.30 3.70 -4.34 7.50 6624 1
## alpha_acs[33] -2.68 3.58 -8.22 3.03 6622 1
## alpha_acs[34] -8.70 3.45 -14.13 -3.17 6123 1
## alpha_acs[35] -6.03 3.56 -11.56 -0.18 6366 1
## alpha_acs[36] -6.46 4.32 -13.46 0.27 7156 1
## alpha_shelf[1] -3.03 1.88 -6.09 -0.25 5866 1
## alpha_shelf[2] 3.68 1.92 0.71 6.60 5938 1
## alpha_shelf[3] -0.09 1.87 -2.93 2.86 5612 1
## alpha_shelf[4] -2.40 1.92 -5.46 0.57 5251 1
## alpha_shelf[5] -1.09 1.94 -3.96 2.09 5320 1
## alpha_shelf[6] 2.70 1.89 -0.14 5.81 5229 1
## beta_sp[1] -3.68 4.99 -11.73 4.19 3941 1
## beta_sp[2] 2.65 4.98 -5.19 10.63 3902 1
## beta_sp[3] 2.80 4.96 -4.82 10.94 3488 1
## beta_sp[4] -11.30 5.03 -19.47 -3.37 4248 1
## beta_sp[5] 6.52 4.99 -1.29 14.46 4143 1
## beta_shade 22.02 5.01 14.11 30.14 3920 1
## beta_sp_shade[1] 5.30 4.61 -2.19 12.54 4287 1
## beta_sp_shade[2] -5.43 4.57 -12.36 2.22 4070 1
## beta_sp_shade[3] -1.81 4.56 -9.18 5.31 4081 1
## beta_sp_shade[4] 5.00 4.72 -2.68 12.50 4382 1
## beta_sp_shade[5] 1.92 4.59 -5.84 8.89 4272 1
## sigma 13.05 0.30 12.59 13.53 9577 1
## sigma_acs 7.57 1.10 5.83 9.22 6869 1
## sigma_shelf 2.86 0.91 1.53 4.18 5378 1compare(m.shade,m.shade.species,m.shade.species.int,m.shade.species.int.acs,m.shade.species.int.acs.shelf)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.int.acs.shelf 8083.3 42.3 0.0 1 56.39 NA
## m.shade.species.int.acs 8124.3 38.7 41.0 0 55.22 12.04
## m.shade.species.int 8334.6 10.6 251.3 0 53.30 30.23
## m.shade.species 8348.1 6.9 264.8 0 53.11 31.20
## m.shade 8458.4 3.2 375.1 0 52.47 38.45coeftab(m.shade.species.int.acs,m.shade.species.int.acs.shelf)## m.shade.species.int.acs m.shade.species.int.acs.shelf
## alpha 41.50 41.28
## alpha_acs[1] -5.40 -6.46
## alpha_acs[2] -10.63 -11.35
## alpha_acs[3] -3.54 -3.32
## alpha_acs[4] 3.15 3.43
## alpha_acs[5] 13.56 13.88
## alpha_acs[6] 2.48 2.73
## alpha_acs[7] 9.86 9.93
## alpha_acs[8] 6.35 6.89
## alpha_acs[9] 3.44 3.63
## alpha_acs[10] 3.86 2.97
## alpha_acs[11] 2.31 2.70
## alpha_acs[12] -9.78 -10.18
## alpha_acs[13] -4.01 -4.09
## alpha_acs[14] -0.69 -0.79
## alpha_acs[15] -3.50 -5.53
## alpha_acs[16] 1.36 1.86
## alpha_acs[17] 4.78 5.00
## alpha_acs[18] -12.07 -12.38
## alpha_acs[19] 15.04 15.11
## alpha_acs[20] -6.18 -5.13
## alpha_acs[21] 0.96 1.60
## alpha_acs[22] 6.18 4.56
## alpha_acs[23] -2.24 -2.21
## alpha_acs[24] -4.28 -4.02
## alpha_acs[25] -5.10 -4.56
## alpha_acs[26] -0.31 -0.19
## alpha_acs[27] 1.83 2.19
## alpha_acs[28] -4.60 -4.73
## alpha_acs[29] 1.14 1.26
## alpha_acs[30] 6.78 6.76
## alpha_acs[31] 10.34 10.76
## alpha_acs[32] 0.51 1.30
## alpha_acs[33] -3.30 -2.68
## alpha_acs[34] -8.61 -8.70
## alpha_acs[35] -5.57 -6.03
## alpha_acs[36] -7.25 -6.46
## beta_sp[1] -4.01 -3.68
## beta_sp[2] 2.54 2.65
## beta_sp[3] 2.65 2.80
## beta_sp[4] -10.69 -11.30
## beta_sp[5] 6.23 6.52
## beta_shade 22.05 22.02
## beta_sp_shade[1] 5.72 5.30
## beta_sp_shade[2] -5.12 -5.43
## beta_sp_shade[3] -1.89 -1.81
## beta_sp_shade[4] 5.11 5.00
## beta_sp_shade[5] 2.39 1.92
## sigma 13.34 13.05
## sigma_acs 7.47 7.57
## alpha_shelf[1] NA -3.03
## alpha_shelf[2] NA 3.68
## alpha_shelf[3] NA -0.09
## alpha_shelf[4] NA -2.4
## alpha_shelf[5] NA -1.09
## alpha_shelf[6] NA 2.7
## sigma_shelf NA 2.86
## nobs 1008 1008plot(coeftab(m.shade.species.int.acs,m.shade.species.int.acs.shelf),pars=c(1,38:50))
including shelf does not change the estimates much, although that model is apparently preferred
reparamterize with no “main” shade effect
m.shade.species.sep.acs.shelf <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_sp_shade[species_id]*shade,
alpha ~ dnorm(53,20),
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(0,10),
beta_sp_shade[species_id] ~ dnorm(0,20),
c(sigma,sigma_acs) ~ dcauchy(0,1),
sigma_shelf ~ dexp(1) # with only 3 shelves per trt might be hard to estimate w. cauchy
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000306 seconds (Sampling)
## 0.00031 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image(file="tomato1212.Rdata")plot(m.shade.species.sep.acs.shelf,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
par(mfrow=c(1,1))
precis(m.shade.species.sep.acs.shelf,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 42.01 4.90 34.24 49.73 2020 1
## alpha_acs[1] -6.49 3.76 -12.47 -0.57 3002 1
## alpha_acs[2] -11.40 3.57 -17.33 -5.94 3545 1
## alpha_acs[3] -3.35 3.43 -8.67 2.17 2703 1
## alpha_acs[4] 3.41 3.43 -1.91 9.13 2713 1
## alpha_acs[5] 13.88 3.49 8.29 19.32 2924 1
## alpha_acs[6] 2.70 3.50 -2.95 8.11 2860 1
## alpha_acs[7] 9.92 3.47 4.12 15.17 3288 1
## alpha_acs[8] 6.91 3.21 1.92 12.16 3134 1
## alpha_acs[9] 3.66 3.27 -1.75 8.70 3300 1
## alpha_acs[10] 2.87 3.94 -3.23 9.20 3897 1
## alpha_acs[11] 2.66 3.36 -2.51 8.11 3133 1
## alpha_acs[12] -10.18 3.42 -15.86 -4.98 3582 1
## alpha_acs[13] -4.14 3.71 -10.25 1.53 3588 1
## alpha_acs[14] -0.82 3.50 -6.40 4.71 3933 1
## alpha_acs[15] -5.58 4.95 -13.28 2.52 5734 1
## alpha_acs[16] 1.78 3.35 -3.72 6.97 2886 1
## alpha_acs[17] 4.99 3.46 -0.34 10.70 3087 1
## alpha_acs[18] -12.46 3.52 -18.18 -7.01 3141 1
## alpha_acs[19] 15.04 3.51 9.30 20.49 2977 1
## alpha_acs[20] -5.13 3.20 -10.40 -0.11 3344 1
## alpha_acs[21] 1.59 3.23 -3.82 6.57 3227 1
## alpha_acs[22] 4.55 4.56 -2.97 11.44 4865 1
## alpha_acs[23] -2.26 3.78 -8.21 3.74 3793 1
## alpha_acs[24] -4.07 3.43 -9.51 1.44 2774 1
## alpha_acs[25] -4.59 3.35 -9.99 0.61 2713 1
## alpha_acs[26] -0.26 3.66 -6.21 5.49 3503 1
## alpha_acs[27] 2.12 3.42 -3.42 7.50 2800 1
## alpha_acs[28] -4.79 3.48 -10.30 0.80 3069 1
## alpha_acs[29] 1.21 3.60 -4.34 7.14 3121 1
## alpha_acs[30] 6.71 3.48 1.15 12.20 3136 1
## alpha_acs[31] 10.75 3.39 5.38 16.27 3440 1
## alpha_acs[32] 1.22 3.69 -4.89 6.88 3554 1
## alpha_acs[33] -2.74 3.64 -8.46 3.13 3490 1
## alpha_acs[34] -8.73 3.48 -14.30 -3.24 3253 1
## alpha_acs[35] -6.04 3.59 -11.70 -0.26 3440 1
## alpha_acs[36] -6.42 4.39 -13.11 0.83 4647 1
## alpha_shelf[1] -2.43 1.96 -5.58 0.56 2130 1
## alpha_shelf[2] 4.35 2.06 1.02 7.47 1925 1
## alpha_shelf[3] 0.55 1.98 -2.62 3.61 1936 1
## alpha_shelf[4] -3.16 2.03 -6.27 0.03 1923 1
## alpha_shelf[5] -1.82 2.02 -4.86 1.41 2074 1
## alpha_shelf[6] 1.97 1.94 -0.85 5.16 2052 1
## beta_sp[1] -3.58 5.01 -11.59 4.46 2411 1
## beta_sp[2] 2.84 5.06 -5.14 10.89 2158 1
## beta_sp[3] 2.90 5.03 -4.88 11.12 2198 1
## beta_sp[4] -11.16 5.20 -19.30 -2.76 2259 1
## beta_sp[5] 6.62 5.14 -1.37 14.84 2173 1
## beta_sp_shade[1] 25.88 3.17 20.71 30.66 1743 1
## beta_sp_shade[2] 14.93 3.11 10.08 19.88 1671 1
## beta_sp_shade[3] 18.64 3.12 13.91 23.65 1691 1
## beta_sp_shade[4] 25.51 3.46 20.00 30.99 2004 1
## beta_sp_shade[5] 22.47 3.14 17.34 27.20 1672 1
## sigma 13.05 0.29 12.55 13.48 8817 1
## sigma_acs 7.57 1.08 5.92 9.25 5285 1
## sigma_shelf 2.99 0.98 1.56 4.38 2973 1compare(m.shade,m.shade.species,m.shade.species.int,m.shade.species.int.acs,m.shade.species.int.acs.shelf,m.shade.species.sep.acs.shelf)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.int.acs.shelf 8083.3 42.3 0.0 0.55 56.39 NA
## m.shade.species.sep.acs.shelf 8083.7 42.5 0.4 0.45 56.40 0.42
## m.shade.species.int.acs 8124.3 38.7 41.0 0.00 55.22 12.04
## m.shade.species.int 8334.6 10.6 251.3 0.00 53.30 30.23
## m.shade.species 8348.1 6.9 264.8 0.00 53.11 31.20
## m.shade 8458.4 3.2 375.1 0.00 52.47 38.45coeftab(m.shade.species.sep.acs.shelf,m.shade.species.int.acs.shelf)## m.shade.species.sep.acs.shelf
## alpha 42.01
## alpha_acs[1] -6.49
## alpha_acs[2] -11.40
## alpha_acs[3] -3.35
## alpha_acs[4] 3.41
## alpha_acs[5] 13.88
## alpha_acs[6] 2.70
## alpha_acs[7] 9.92
## alpha_acs[8] 6.91
## alpha_acs[9] 3.66
## alpha_acs[10] 2.87
## alpha_acs[11] 2.66
## alpha_acs[12] -10.18
## alpha_acs[13] -4.14
## alpha_acs[14] -0.82
## alpha_acs[15] -5.58
## alpha_acs[16] 1.78
## alpha_acs[17] 4.99
## alpha_acs[18] -12.46
## alpha_acs[19] 15.04
## alpha_acs[20] -5.13
## alpha_acs[21] 1.59
## alpha_acs[22] 4.55
## alpha_acs[23] -2.26
## alpha_acs[24] -4.07
## alpha_acs[25] -4.59
## alpha_acs[26] -0.26
## alpha_acs[27] 2.12
## alpha_acs[28] -4.79
## alpha_acs[29] 1.21
## alpha_acs[30] 6.71
## alpha_acs[31] 10.75
## alpha_acs[32] 1.22
## alpha_acs[33] -2.74
## alpha_acs[34] -8.73
## alpha_acs[35] -6.04
## alpha_acs[36] -6.42
## alpha_shelf[1] -2.43
## alpha_shelf[2] 4.35
## alpha_shelf[3] 0.55
## alpha_shelf[4] -3.16
## alpha_shelf[5] -1.82
## alpha_shelf[6] 1.97
## beta_sp[1] -3.58
## beta_sp[2] 2.84
## beta_sp[3] 2.9
## beta_sp[4] -11.16
## beta_sp[5] 6.62
## beta_sp_shade[1] 25.88
## beta_sp_shade[2] 14.93
## beta_sp_shade[3] 18.64
## beta_sp_shade[4] 25.51
## beta_sp_shade[5] 22.47
## sigma 13.05
## sigma_acs 7.57
## sigma_shelf 2.99
## beta_shade NA
## nobs 1008
## m.shade.species.int.acs.shelf
## alpha 41.28
## alpha_acs[1] -6.46
## alpha_acs[2] -11.35
## alpha_acs[3] -3.32
## alpha_acs[4] 3.43
## alpha_acs[5] 13.88
## alpha_acs[6] 2.73
## alpha_acs[7] 9.93
## alpha_acs[8] 6.89
## alpha_acs[9] 3.63
## alpha_acs[10] 2.97
## alpha_acs[11] 2.70
## alpha_acs[12] -10.18
## alpha_acs[13] -4.09
## alpha_acs[14] -0.79
## alpha_acs[15] -5.53
## alpha_acs[16] 1.86
## alpha_acs[17] 5.00
## alpha_acs[18] -12.38
## alpha_acs[19] 15.11
## alpha_acs[20] -5.13
## alpha_acs[21] 1.60
## alpha_acs[22] 4.56
## alpha_acs[23] -2.21
## alpha_acs[24] -4.02
## alpha_acs[25] -4.56
## alpha_acs[26] -0.19
## alpha_acs[27] 2.19
## alpha_acs[28] -4.73
## alpha_acs[29] 1.26
## alpha_acs[30] 6.76
## alpha_acs[31] 10.76
## alpha_acs[32] 1.30
## alpha_acs[33] -2.68
## alpha_acs[34] -8.70
## alpha_acs[35] -6.03
## alpha_acs[36] -6.46
## alpha_shelf[1] -3.03
## alpha_shelf[2] 3.68
## alpha_shelf[3] -0.09
## alpha_shelf[4] -2.40
## alpha_shelf[5] -1.09
## alpha_shelf[6] 2.70
## beta_sp[1] -3.68
## beta_sp[2] 2.65
## beta_sp[3] 2.8
## beta_sp[4] -11.30
## beta_sp[5] 6.52
## beta_sp_shade[1] 5.30
## beta_sp_shade[2] -5.43
## beta_sp_shade[3] -1.81
## beta_sp_shade[4] 5.00
## beta_sp_shade[5] 1.92
## sigma 13.05
## sigma_acs 7.57
## sigma_shelf 2.86
## beta_shade 22.02
## nobs 1008plot(coeftab(m.shade.species.sep.acs.shelf))
This parameterization leads to the same WAIC, so it is an equivalent model. Interestingly though, the SD on the shade coefficients is smaller, presumably because one fewer parameter is being estimated.
Can we do no main alpha?
m.shade.species.sep.acs.shelf.no.alpha <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <-
alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_sp_shade[species_id]*shade,
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(42,20),
beta_sp_shade[species_id] ~ dnorm(0,20),
c(sigma,sigma_acs) ~ dcauchy(0,1),
sigma_shelf ~ dexp(1) # with only 3 shelves per trt might be hard to estimate w. cauchy
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 6e-06 seconds (Warm-up)
## 0.000312 seconds (Sampling)
## 0.000318 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image(file="tomato1212.Rdata")plot(m.shade.species.sep.acs.shelf.no.alpha,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
par(mfrow=c(1,1))
precis(m.shade.species.sep.acs.shelf.no.alpha,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_acs[1] -6.79 3.75 -12.75 -0.89 3134 1
## alpha_acs[2] -11.67 3.64 -17.38 -5.74 3279 1
## alpha_acs[3] -3.60 3.48 -9.29 1.83 2720 1
## alpha_acs[4] 3.13 3.47 -2.59 8.43 2819 1
## alpha_acs[5] 13.62 3.47 7.92 18.95 2685 1
## alpha_acs[6] 2.41 3.48 -3.06 8.08 2735 1
## alpha_acs[7] 9.59 3.52 4.14 15.31 2834 1
## alpha_acs[8] 6.75 3.34 1.18 11.88 3047 1
## alpha_acs[9] 3.53 3.38 -1.99 8.78 3022 1
## alpha_acs[10] 3.75 3.99 -2.58 10.09 3719 1
## alpha_acs[11] 2.34 3.46 -3.11 7.84 2778 1
## alpha_acs[12] -10.30 3.52 -15.74 -4.51 3340 1
## alpha_acs[13] -3.29 3.76 -9.11 2.79 3303 1
## alpha_acs[14] -0.95 3.63 -6.66 4.94 3578 1
## alpha_acs[15] -5.61 4.93 -13.48 2.27 6084 1
## alpha_acs[16] 1.48 3.39 -3.70 7.11 2573 1
## alpha_acs[17] 5.14 3.48 -0.59 10.54 2873 1
## alpha_acs[18] -12.26 3.62 -17.81 -6.37 3211 1
## alpha_acs[19] 14.69 3.57 9.26 20.53 2797 1
## alpha_acs[20] -5.27 3.34 -10.70 -0.03 2946 1
## alpha_acs[21] 1.46 3.38 -3.93 6.84 3113 1
## alpha_acs[22] 5.21 4.57 -2.08 12.54 4562 1
## alpha_acs[23] -1.44 3.82 -7.39 4.78 3462 1
## alpha_acs[24] -4.32 3.45 -9.58 1.43 2600 1
## alpha_acs[25] -4.84 3.42 -10.35 0.59 2649 1
## alpha_acs[26] -0.09 3.69 -5.94 5.81 3262 1
## alpha_acs[27] 2.28 3.44 -3.45 7.49 2836 1
## alpha_acs[28] -4.63 3.52 -10.19 1.04 2968 1
## alpha_acs[29] 1.38 3.60 -4.47 6.97 3066 1
## alpha_acs[30] 6.85 3.47 1.28 12.28 2882 1
## alpha_acs[31] 10.66 3.49 5.08 16.15 3252 1
## alpha_acs[32] 2.11 3.71 -3.79 8.02 3079 1
## alpha_acs[33] -1.90 3.66 -7.62 4.01 3242 1
## alpha_acs[34] -9.04 3.51 -14.89 -3.72 2867 1
## alpha_acs[35] -6.39 3.60 -12.04 -0.66 3076 1
## alpha_acs[36] -5.70 4.38 -12.58 1.37 4809 1
## alpha_shelf[1] -2.38 1.96 -5.57 0.60 2963 1
## alpha_shelf[2] 4.41 2.06 1.35 7.73 2804 1
## alpha_shelf[3] 0.60 1.97 -2.36 3.78 2811 1
## alpha_shelf[4] -3.05 1.96 -5.99 0.12 1925 1
## alpha_shelf[5] -1.68 1.95 -4.77 1.37 2061 1
## alpha_shelf[6] 2.10 1.88 -0.89 5.00 2010 1
## beta_sp[1] 38.09 3.49 32.76 43.85 1965 1
## beta_sp[2] 45.03 3.55 39.32 50.61 1930 1
## beta_sp[3] 44.91 3.38 39.77 50.58 1980 1
## beta_sp[4] 29.62 3.74 23.65 35.58 2207 1
## beta_sp[5] 48.90 3.55 43.22 54.42 1953 1
## beta_sp_shade[1] 26.03 3.08 21.38 31.05 2237 1
## beta_sp_shade[2] 14.96 3.04 10.49 20.03 2169 1
## beta_sp_shade[3] 18.69 3.02 14.13 23.70 2219 1
## beta_sp_shade[4] 25.85 3.34 20.73 31.29 2538 1
## beta_sp_shade[5] 22.44 3.05 17.62 27.21 2178 1
## sigma 13.05 0.30 12.58 13.53 10166 1
## sigma_acs 7.57 1.09 5.81 9.16 5217 1
## sigma_shelf 2.98 0.96 1.57 4.34 3519 1compare(m.shade,m.shade.species,m.shade.species.int,m.shade.species.int.acs,m.shade.species.int.acs.shelf,m.shade.species.sep.acs.shelf,m.shade.species.sep.acs.shelf.no.alpha)## WAIC pWAIC dWAIC weight SE
## m.shade.species.sep.acs.shelf.no.alpha 8083.3 42.3 0.0 0.36 56.42
## m.shade.species.int.acs.shelf 8083.3 42.3 0.0 0.35 56.39
## m.shade.species.sep.acs.shelf 8083.7 42.5 0.5 0.28 56.40
## m.shade.species.int.acs 8124.3 38.7 41.0 0.00 55.22
## m.shade.species.int 8334.6 10.6 251.4 0.00 53.30
## m.shade.species 8348.1 6.9 264.9 0.00 53.11
## m.shade 8458.4 3.2 375.1 0.00 52.47
## dSE
## m.shade.species.sep.acs.shelf.no.alpha NA
## m.shade.species.int.acs.shelf 0.54
## m.shade.species.sep.acs.shelf 0.32
## m.shade.species.int.acs 12.17
## m.shade.species.int 30.19
## m.shade.species 31.19
## m.shade 38.46coeftab(m.shade.species.sep.acs.shelf,m.shade.species.int.acs.shelf)## m.shade.species.sep.acs.shelf
## alpha 42.01
## alpha_acs[1] -6.49
## alpha_acs[2] -11.40
## alpha_acs[3] -3.35
## alpha_acs[4] 3.41
## alpha_acs[5] 13.88
## alpha_acs[6] 2.70
## alpha_acs[7] 9.92
## alpha_acs[8] 6.91
## alpha_acs[9] 3.66
## alpha_acs[10] 2.87
## alpha_acs[11] 2.66
## alpha_acs[12] -10.18
## alpha_acs[13] -4.14
## alpha_acs[14] -0.82
## alpha_acs[15] -5.58
## alpha_acs[16] 1.78
## alpha_acs[17] 4.99
## alpha_acs[18] -12.46
## alpha_acs[19] 15.04
## alpha_acs[20] -5.13
## alpha_acs[21] 1.59
## alpha_acs[22] 4.55
## alpha_acs[23] -2.26
## alpha_acs[24] -4.07
## alpha_acs[25] -4.59
## alpha_acs[26] -0.26
## alpha_acs[27] 2.12
## alpha_acs[28] -4.79
## alpha_acs[29] 1.21
## alpha_acs[30] 6.71
## alpha_acs[31] 10.75
## alpha_acs[32] 1.22
## alpha_acs[33] -2.74
## alpha_acs[34] -8.73
## alpha_acs[35] -6.04
## alpha_acs[36] -6.42
## alpha_shelf[1] -2.43
## alpha_shelf[2] 4.35
## alpha_shelf[3] 0.55
## alpha_shelf[4] -3.16
## alpha_shelf[5] -1.82
## alpha_shelf[6] 1.97
## beta_sp[1] -3.58
## beta_sp[2] 2.84
## beta_sp[3] 2.9
## beta_sp[4] -11.16
## beta_sp[5] 6.62
## beta_sp_shade[1] 25.88
## beta_sp_shade[2] 14.93
## beta_sp_shade[3] 18.64
## beta_sp_shade[4] 25.51
## beta_sp_shade[5] 22.47
## sigma 13.05
## sigma_acs 7.57
## sigma_shelf 2.99
## beta_shade NA
## nobs 1008
## m.shade.species.int.acs.shelf
## alpha 41.28
## alpha_acs[1] -6.46
## alpha_acs[2] -11.35
## alpha_acs[3] -3.32
## alpha_acs[4] 3.43
## alpha_acs[5] 13.88
## alpha_acs[6] 2.73
## alpha_acs[7] 9.93
## alpha_acs[8] 6.89
## alpha_acs[9] 3.63
## alpha_acs[10] 2.97
## alpha_acs[11] 2.70
## alpha_acs[12] -10.18
## alpha_acs[13] -4.09
## alpha_acs[14] -0.79
## alpha_acs[15] -5.53
## alpha_acs[16] 1.86
## alpha_acs[17] 5.00
## alpha_acs[18] -12.38
## alpha_acs[19] 15.11
## alpha_acs[20] -5.13
## alpha_acs[21] 1.60
## alpha_acs[22] 4.56
## alpha_acs[23] -2.21
## alpha_acs[24] -4.02
## alpha_acs[25] -4.56
## alpha_acs[26] -0.19
## alpha_acs[27] 2.19
## alpha_acs[28] -4.73
## alpha_acs[29] 1.26
## alpha_acs[30] 6.76
## alpha_acs[31] 10.76
## alpha_acs[32] 1.30
## alpha_acs[33] -2.68
## alpha_acs[34] -8.70
## alpha_acs[35] -6.03
## alpha_acs[36] -6.46
## alpha_shelf[1] -3.03
## alpha_shelf[2] 3.68
## alpha_shelf[3] -0.09
## alpha_shelf[4] -2.40
## alpha_shelf[5] -1.09
## alpha_shelf[6] 2.70
## beta_sp[1] -3.68
## beta_sp[2] 2.65
## beta_sp[3] 2.8
## beta_sp[4] -11.30
## beta_sp[5] 6.52
## beta_sp_shade[1] 5.30
## beta_sp_shade[2] -5.43
## beta_sp_shade[3] -1.81
## beta_sp_shade[4] 5.00
## beta_sp_shade[5] 1.92
## sigma 13.05
## sigma_acs 7.57
## sigma_shelf 2.86
## beta_shade 22.02
## nobs 1008plot(coeftab(m.shade.species.sep.acs.shelf.no.alpha))
So long as an appropriate prior for the species intercepts is used this is OK. The other reasonable paramaterization would be to have the alpha befor the first species.
Consider different slopes for accessions
m.shade.species.sep.acs.int.shelf.no.alpha <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_acs_shade[acs_id]*shade,
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(42,10),
beta_acs_shade[acs_id] ~ dnorm(20,20),
c(sigma,sigma_acs,sigma_shelf) ~ dcauchy(0,1)
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.000412 seconds (Sampling)
## 0.000415 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image("tomato1212.Rdata")plot(m.shade.species.sep.acs.int.shelf.no.alpha,ask=FALSE)## Waiting to draw page 2 of 6
## Waiting to draw page 3 of 6
## Waiting to draw page 4 of 6
## Waiting to draw page 5 of 6
## Waiting to draw page 6 of 6
par(mfrow=c(1,1))
precis(m.shade.species.sep.acs.int.shelf.no.alpha,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_acs[1] -1.88 3.27 -7.08 3.26 12000 1
## alpha_acs[2] -4.05 3.29 -9.52 1.06 12000 1
## alpha_acs[3] -4.01 3.02 -8.77 0.88 12000 1
## alpha_acs[4] 2.83 3.05 -1.90 7.77 12000 1
## alpha_acs[5] 8.06 3.17 2.85 12.93 12000 1
## alpha_acs[6] -0.15 3.01 -5.00 4.57 12000 1
## alpha_acs[7] 5.86 3.18 0.88 10.96 12000 1
## alpha_acs[8] 0.12 2.86 -4.56 4.54 12000 1
## alpha_acs[9] 2.53 2.95 -2.13 7.23 12000 1
## alpha_acs[10] -0.09 3.42 -5.69 5.19 12000 1
## alpha_acs[11] 3.03 2.96 -1.82 7.65 12000 1
## alpha_acs[12] -4.98 3.12 -9.93 0.06 12000 1
## alpha_acs[13] 0.25 3.32 -5.26 5.43 12000 1
## alpha_acs[14] 2.08 3.11 -3.05 6.82 12000 1
## alpha_acs[15] -1.06 3.95 -7.59 4.98 12000 1
## alpha_acs[16] 2.69 2.88 -2.07 7.10 12000 1
## alpha_acs[17] 4.76 3.03 0.22 9.78 12000 1
## alpha_acs[18] -7.57 3.21 -12.78 -2.63 12000 1
## alpha_acs[19] 5.15 3.16 0.24 10.35 12000 1
## alpha_acs[20] -3.38 2.89 -8.12 1.01 12000 1
## alpha_acs[21] -0.07 2.85 -4.51 4.55 12000 1
## alpha_acs[22] 0.55 3.95 -5.38 7.24 12000 1
## alpha_acs[23] -1.77 3.42 -7.32 3.57 12000 1
## alpha_acs[24] -2.48 2.96 -7.32 2.11 12000 1
## alpha_acs[25] -1.83 2.91 -6.49 2.76 12000 1
## alpha_acs[26] -2.09 3.29 -7.23 3.31 12000 1
## alpha_acs[27] 1.02 2.92 -3.59 5.71 12000 1
## alpha_acs[28] -0.43 3.12 -5.41 4.51 12000 1
## alpha_acs[29] -0.57 3.12 -5.44 4.50 12000 1
## alpha_acs[30] 3.83 3.12 -1.21 8.75 12000 1
## alpha_acs[31] 5.51 3.17 0.45 10.55 12000 1
## alpha_acs[32] -0.14 3.14 -5.00 4.96 12000 1
## alpha_acs[33] 0.40 3.25 -4.93 5.40 12000 1
## alpha_acs[34] -6.86 3.14 -11.89 -1.88 12000 1
## alpha_acs[35] -4.02 3.13 -8.91 1.03 12000 1
## alpha_acs[36] -2.10 3.88 -8.21 4.05 12000 1
## alpha_shelf[1] -3.10 1.91 -6.14 -0.15 6532 1
## alpha_shelf[2] 4.20 1.92 1.15 7.16 6021 1
## alpha_shelf[3] 0.41 1.88 -2.56 3.33 6135 1
## alpha_shelf[4] -2.56 1.78 -5.36 0.25 5320 1
## alpha_shelf[5] -1.56 1.81 -4.18 1.50 5755 1
## alpha_shelf[6] 2.30 1.73 -0.28 5.16 5317 1
## beta_sp[1] 37.82 2.60 33.70 41.95 6643 1
## beta_sp[2] 44.54 2.56 40.49 48.58 7021 1
## beta_sp[3] 45.15 2.53 40.98 48.96 5905 1
## beta_sp[4] 29.80 2.78 25.65 34.53 7823 1
## beta_sp[5] 48.55 2.58 44.44 52.66 6719 1
## beta_acs_shade[1] 5.52 5.34 -3.22 13.76 12000 1
## beta_acs_shade[2] 9.18 4.93 1.15 16.87 8680 1
## beta_acs_shade[3] 17.79 4.39 11.07 25.17 12000 1
## beta_acs_shade[4] 15.83 4.45 8.40 22.67 12000 1
## beta_acs_shade[5] 25.38 4.51 17.87 32.31 8931 1
## beta_acs_shade[6] 21.93 4.58 14.74 29.19 12000 1
## beta_acs_shade[7] 29.12 4.55 22.18 36.64 12000 1
## beta_acs_shade[8] 33.92 4.33 26.84 40.65 12000 1
## beta_acs_shade[9] 19.81 4.40 12.96 26.95 12000 1
## beta_acs_shade[10] 34.62 5.65 25.42 43.49 12000 1
## beta_acs_shade[11] 21.04 4.27 13.82 27.46 12000 1
## beta_acs_shade[12] 7.13 4.89 -0.61 14.91 12000 1
## beta_acs_shade[13] 18.76 4.84 11.32 26.87 12000 1
## beta_acs_shade[14] 9.22 5.10 1.01 17.28 12000 1
## beta_acs_shade[15] 1.91 8.92 -11.61 16.61 12000 1
## beta_acs_shade[16] 20.01 4.14 13.71 26.88 9347 1
## beta_acs_shade[17] 26.18 4.49 19.47 33.80 12000 1
## beta_acs_shade[18] 18.01 4.84 10.34 25.81 12000 1
## beta_acs_shade[19] 41.62 4.64 34.17 48.97 12000 1
## beta_acs_shade[20] 14.98 4.30 8.33 22.08 12000 1
## beta_acs_shade[21] 22.02 4.37 15.24 29.14 12000 1
## beta_acs_shade[22] 35.36 6.95 24.57 46.50 12000 1
## beta_acs_shade[23] 27.11 5.14 19.04 35.33 12000 1
## beta_acs_shade[24] 12.70 4.32 5.95 19.66 12000 1
## beta_acs_shade[25] 10.05 4.16 3.54 16.82 8890 1
## beta_acs_shade[26] 31.63 5.12 23.35 39.70 12000 1
## beta_acs_shade[27] 29.20 4.27 22.37 35.99 12000 1
## beta_acs_shade[28] 17.83 4.64 10.36 25.03 12000 1
## beta_acs_shade[29] 31.18 4.85 23.25 38.76 12000 1
## beta_acs_shade[30] 31.55 4.43 24.31 38.43 12000 1
## beta_acs_shade[31] 27.45 4.63 20.25 35.05 12000 1
## beta_acs_shade[32] 31.17 4.95 23.48 39.18 12000 1
## beta_acs_shade[33] 21.07 4.70 13.45 28.37 12000 1
## beta_acs_shade[34] 20.51 4.64 12.99 27.86 12000 1
## beta_acs_shade[35] 19.18 5.01 11.42 27.42 12000 1
## beta_acs_shade[36] 19.26 6.52 8.91 29.80 12000 1
## sigma 12.73 0.29 12.28 13.21 12000 1
## sigma_acs 4.74 0.91 3.26 6.12 4913 1
## sigma_shelf 3.32 1.32 1.57 4.99 6746 1compare(m.shade.species.sep.acs.int.shelf.no.alpha,m.shade.species.sep.acs.shelf.no.alpha)## WAIC pWAIC dWAIC weight SE
## m.shade.species.sep.acs.int.shelf.no.alpha 8062.2 67.3 0.0 1 58.73
## m.shade.species.sep.acs.shelf.no.alpha 8083.3 42.3 21.1 0 56.42
## dSE
## m.shade.species.sep.acs.int.shelf.no.alpha NA
## m.shade.species.sep.acs.shelf.no.alpha 20.36coeftab(m.shade.species.sep.acs.int.shelf.no.alpha,m.shade.species.sep.acs.shelf.no.alpha)## m.shade.species.sep.acs.int.shelf.no.alpha
## alpha_acs[1] -1.88
## alpha_acs[2] -4.05
## alpha_acs[3] -4.01
## alpha_acs[4] 2.83
## alpha_acs[5] 8.06
## alpha_acs[6] -0.15
## alpha_acs[7] 5.86
## alpha_acs[8] 0.12
## alpha_acs[9] 2.53
## alpha_acs[10] -0.09
## alpha_acs[11] 3.03
## alpha_acs[12] -4.98
## alpha_acs[13] 0.25
## alpha_acs[14] 2.08
## alpha_acs[15] -1.06
## alpha_acs[16] 2.69
## alpha_acs[17] 4.76
## alpha_acs[18] -7.57
## alpha_acs[19] 5.15
## alpha_acs[20] -3.38
## alpha_acs[21] -0.07
## alpha_acs[22] 0.55
## alpha_acs[23] -1.77
## alpha_acs[24] -2.48
## alpha_acs[25] -1.83
## alpha_acs[26] -2.09
## alpha_acs[27] 1.02
## alpha_acs[28] -0.43
## alpha_acs[29] -0.57
## alpha_acs[30] 3.83
## alpha_acs[31] 5.51
## alpha_acs[32] -0.14
## alpha_acs[33] 0.4
## alpha_acs[34] -6.86
## alpha_acs[35] -4.02
## alpha_acs[36] -2.1
## alpha_shelf[1] -3.10
## alpha_shelf[2] 4.20
## alpha_shelf[3] 0.41
## alpha_shelf[4] -2.56
## alpha_shelf[5] -1.56
## alpha_shelf[6] 2.3
## beta_sp[1] 37.82
## beta_sp[2] 44.54
## beta_sp[3] 45.15
## beta_sp[4] 29.80
## beta_sp[5] 48.55
## beta_acs_shade[1] 5.52
## beta_acs_shade[2] 9.18
## beta_acs_shade[3] 17.79
## beta_acs_shade[4] 15.83
## beta_acs_shade[5] 25.38
## beta_acs_shade[6] 21.93
## beta_acs_shade[7] 29.12
## beta_acs_shade[8] 33.92
## beta_acs_shade[9] 19.81
## beta_acs_shade[10] 34.62
## beta_acs_shade[11] 21.04
## beta_acs_shade[12] 7.13
## beta_acs_shade[13] 18.76
## beta_acs_shade[14] 9.22
## beta_acs_shade[15] 1.91
## beta_acs_shade[16] 20.01
## beta_acs_shade[17] 26.18
## beta_acs_shade[18] 18.01
## beta_acs_shade[19] 41.62
## beta_acs_shade[20] 14.98
## beta_acs_shade[21] 22.02
## beta_acs_shade[22] 35.36
## beta_acs_shade[23] 27.11
## beta_acs_shade[24] 12.7
## beta_acs_shade[25] 10.05
## beta_acs_shade[26] 31.63
## beta_acs_shade[27] 29.2
## beta_acs_shade[28] 17.83
## beta_acs_shade[29] 31.18
## beta_acs_shade[30] 31.55
## beta_acs_shade[31] 27.45
## beta_acs_shade[32] 31.17
## beta_acs_shade[33] 21.07
## beta_acs_shade[34] 20.51
## beta_acs_shade[35] 19.18
## beta_acs_shade[36] 19.26
## sigma 12.73
## sigma_acs 4.74
## sigma_shelf 3.32
## beta_sp_shade[1] NA
## beta_sp_shade[2] NA
## beta_sp_shade[3] NA
## beta_sp_shade[4] NA
## beta_sp_shade[5] NA
## nobs 1008
## m.shade.species.sep.acs.shelf.no.alpha
## alpha_acs[1] -6.79
## alpha_acs[2] -11.67
## alpha_acs[3] -3.60
## alpha_acs[4] 3.13
## alpha_acs[5] 13.62
## alpha_acs[6] 2.41
## alpha_acs[7] 9.59
## alpha_acs[8] 6.75
## alpha_acs[9] 3.53
## alpha_acs[10] 3.75
## alpha_acs[11] 2.34
## alpha_acs[12] -10.30
## alpha_acs[13] -3.29
## alpha_acs[14] -0.95
## alpha_acs[15] -5.61
## alpha_acs[16] 1.48
## alpha_acs[17] 5.14
## alpha_acs[18] -12.26
## alpha_acs[19] 14.69
## alpha_acs[20] -5.27
## alpha_acs[21] 1.46
## alpha_acs[22] 5.21
## alpha_acs[23] -1.44
## alpha_acs[24] -4.32
## alpha_acs[25] -4.84
## alpha_acs[26] -0.09
## alpha_acs[27] 2.28
## alpha_acs[28] -4.63
## alpha_acs[29] 1.38
## alpha_acs[30] 6.85
## alpha_acs[31] 10.66
## alpha_acs[32] 2.11
## alpha_acs[33] -1.9
## alpha_acs[34] -9.04
## alpha_acs[35] -6.39
## alpha_acs[36] -5.7
## alpha_shelf[1] -2.38
## alpha_shelf[2] 4.41
## alpha_shelf[3] 0.60
## alpha_shelf[4] -3.05
## alpha_shelf[5] -1.68
## alpha_shelf[6] 2.1
## beta_sp[1] 38.09
## beta_sp[2] 45.03
## beta_sp[3] 44.91
## beta_sp[4] 29.62
## beta_sp[5] 48.90
## beta_acs_shade[1] NA
## beta_acs_shade[2] NA
## beta_acs_shade[3] NA
## beta_acs_shade[4] NA
## beta_acs_shade[5] NA
## beta_acs_shade[6] NA
## beta_acs_shade[7] NA
## beta_acs_shade[8] NA
## beta_acs_shade[9] NA
## beta_acs_shade[10] NA
## beta_acs_shade[11] NA
## beta_acs_shade[12] NA
## beta_acs_shade[13] NA
## beta_acs_shade[14] NA
## beta_acs_shade[15] NA
## beta_acs_shade[16] NA
## beta_acs_shade[17] NA
## beta_acs_shade[18] NA
## beta_acs_shade[19] NA
## beta_acs_shade[20] NA
## beta_acs_shade[21] NA
## beta_acs_shade[22] NA
## beta_acs_shade[23] NA
## beta_acs_shade[24] NA
## beta_acs_shade[25] NA
## beta_acs_shade[26] NA
## beta_acs_shade[27] NA
## beta_acs_shade[28] NA
## beta_acs_shade[29] NA
## beta_acs_shade[30] NA
## beta_acs_shade[31] NA
## beta_acs_shade[32] NA
## beta_acs_shade[33] NA
## beta_acs_shade[34] NA
## beta_acs_shade[35] NA
## beta_acs_shade[36] NA
## sigma 13.05
## sigma_acs 7.57
## sigma_shelf 2.98
## beta_sp_shade[1] 26.03
## beta_sp_shade[2] 14.96
## beta_sp_shade[3] 18.69
## beta_sp_shade[4] 25.85
## beta_sp_shade[5] 22.44
## nobs 1008plot(coeftab(m.shade.species.sep.acs.int.shelf.no.alpha,m.shade.species.sep.acs.shelf.no.alpha),pars=c(1,44:56))
drop the species shade interaction
m.shade.species.acs.int.shelf <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_shade*shade +
beta_acs_shade[acs_id]*shade,
alpha ~ dnorm(53,20),
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
beta_acs_shade[acs_id] ~ dnorm(0,10),
c(sigma,sigma_acs,sigma_shelf) ~ dcauchy(0,1)
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 6e-06 seconds (Warm-up)
## 0.000404 seconds (Sampling)
## 0.00041 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image("tomato1212.Rdata")plot(m.shade.species.acs.int.shelf,ask=FALSE)## Waiting to draw page 2 of 6
## Waiting to draw page 3 of 6
## Waiting to draw page 4 of 6
## Waiting to draw page 5 of 6
## Waiting to draw page 6 of 6
par(mfrow=c(1,1))
precis(m.shade.species.acs.int.shelf,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.50 4.95 33.43 49.17 4629 1
## alpha_acs[1] -2.65 3.32 -7.92 2.70 12000 1
## alpha_acs[2] -4.97 3.28 -10.14 0.28 12000 1
## alpha_acs[3] -4.00 3.12 -8.92 0.99 12000 1
## alpha_acs[4] 3.02 3.13 -1.97 7.90 12000 1
## alpha_acs[5] 8.86 3.18 3.52 13.65 12000 1
## alpha_acs[6] 0.18 3.10 -4.69 5.09 12000 1
## alpha_acs[7] 6.43 3.19 1.13 11.28 12000 1
## alpha_acs[8] 0.88 2.93 -3.79 5.46 12000 1
## alpha_acs[9] 2.83 3.01 -1.78 7.80 12000 1
## alpha_acs[10] 0.39 3.49 -5.07 5.98 12000 1
## alpha_acs[11] 3.07 3.02 -1.69 8.04 12000 1
## alpha_acs[12] -5.60 3.16 -10.64 -0.62 12000 1
## alpha_acs[13] -0.32 3.44 -5.78 5.17 12000 1
## alpha_acs[14] 1.88 3.20 -3.07 7.13 12000 1
## alpha_acs[15] -2.30 4.10 -8.56 4.30 12000 1
## alpha_acs[16] 2.69 2.99 -2.24 7.26 12000 1
## alpha_acs[17] 4.89 3.05 0.21 9.92 12000 1
## alpha_acs[18] -8.27 3.21 -13.04 -2.76 12000 1
## alpha_acs[19] 6.31 3.24 1.33 11.72 12000 1
## alpha_acs[20] -3.45 2.98 -8.18 1.34 12000 1
## alpha_acs[21] 0.24 2.91 -4.53 4.75 12000 1
## alpha_acs[22] 1.56 4.02 -4.76 7.96 12000 1
## alpha_acs[23] -1.91 3.51 -7.41 3.63 12000 1
## alpha_acs[24] -2.62 3.04 -7.48 2.21 12000 1
## alpha_acs[25] -2.02 2.98 -6.91 2.58 12000 1
## alpha_acs[26] -1.83 3.35 -7.12 3.59 12000 1
## alpha_acs[27] 1.13 2.98 -3.69 5.84 12000 1
## alpha_acs[28] -0.90 3.13 -5.82 4.15 12000 1
## alpha_acs[29] -0.36 3.23 -5.32 4.90 12000 1
## alpha_acs[30] 4.25 3.19 -0.75 9.44 12000 1
## alpha_acs[31] 6.29 3.27 1.19 11.51 12000 1
## alpha_acs[32] -0.07 3.25 -5.01 5.35 12000 1
## alpha_acs[33] 0.04 3.30 -5.05 5.43 12000 1
## alpha_acs[34] -7.22 3.18 -12.35 -2.27 12000 1
## alpha_acs[35] -4.39 3.19 -9.50 0.57 12000 1
## alpha_acs[36] -2.68 3.97 -8.92 3.69 12000 1
## alpha_shelf[1] -3.33 2.34 -6.92 0.08 6206 1
## alpha_shelf[2] 4.00 2.38 0.27 7.46 6147 1
## alpha_shelf[3] 0.13 2.34 -3.46 3.57 6271 1
## alpha_shelf[4] -2.40 2.38 -5.97 1.15 5171 1
## alpha_shelf[5] -1.34 2.35 -4.91 2.20 5352 1
## alpha_shelf[6] 2.54 2.33 -1.02 6.05 5418 1
## beta_sp[1] -3.61 4.80 -11.19 4.04 5274 1
## beta_sp[2] 2.37 4.81 -5.23 10.18 5433 1
## beta_sp[3] 3.01 4.77 -4.67 10.50 5439 1
## beta_sp[4] -11.41 4.87 -19.16 -3.72 5825 1
## beta_sp[5] 6.84 4.80 -0.95 14.49 5381 1
## beta_shade 22.08 3.61 16.73 27.73 4417 1
## beta_acs_shade[1] -13.40 4.85 -21.04 -5.57 12000 1
## beta_acs_shade[2] -10.71 4.51 -18.12 -3.71 12000 1
## beta_acs_shade[3] -3.19 4.06 -9.81 3.06 12000 1
## beta_acs_shade[4] -5.14 4.14 -11.73 1.48 12000 1
## beta_acs_shade[5] 3.22 4.13 -3.37 9.77 12000 1
## beta_acs_shade[6] 0.45 4.20 -6.38 6.90 12000 1
## beta_acs_shade[7] 6.54 4.17 0.13 13.38 12000 1
## beta_acs_shade[8] 11.25 4.08 5.21 18.24 12000 1
## beta_acs_shade[9] -1.50 4.07 -8.03 5.03 12000 1
## beta_acs_shade[10] 10.56 5.13 1.99 18.44 12000 1
## beta_acs_shade[11] -0.61 4.03 -7.27 5.60 12000 1
## beta_acs_shade[12] -12.25 4.53 -19.34 -4.94 12000 1
## beta_acs_shade[13] -2.57 4.47 -9.48 4.71 12000 1
## beta_acs_shade[14] -10.53 4.73 -17.90 -2.78 12000 1
## beta_acs_shade[15] -11.97 7.18 -23.52 -0.46 12000 1
## beta_acs_shade[16] -1.49 3.92 -7.62 4.79 12000 1
## beta_acs_shade[17] 3.94 4.22 -2.74 10.64 12000 1
## beta_acs_shade[18] -2.93 4.44 -10.30 3.83 12000 1
## beta_acs_shade[19] 17.64 4.27 11.04 24.70 12000 1
## beta_acs_shade[20] -5.75 4.04 -12.21 0.69 12000 1
## beta_acs_shade[21] 0.53 4.03 -5.87 6.98 12000 1
## beta_acs_shade[22] 10.26 5.90 0.82 19.52 12000 1
## beta_acs_shade[23] 4.65 4.76 -2.88 12.43 12000 1
## beta_acs_shade[24] -7.81 4.03 -14.56 -1.69 12000 1
## beta_acs_shade[25] -10.29 3.93 -16.37 -3.85 12000 1
## beta_acs_shade[26] 8.61 4.55 1.56 15.98 12000 1
## beta_acs_shade[27] 6.76 4.05 0.44 13.24 12000 1
## beta_acs_shade[28] -3.36 4.25 -10.21 3.37 12000 1
## beta_acs_shade[29] 8.40 4.46 1.30 15.58 12000 1
## beta_acs_shade[30] 8.65 4.10 2.22 15.29 12000 1
## beta_acs_shade[31] 5.16 4.40 -1.94 12.14 12000 1
## beta_acs_shade[32] 8.15 4.54 1.15 15.56 12000 1
## beta_acs_shade[33] -0.57 4.29 -7.48 6.23 12000 1
## beta_acs_shade[34] -0.70 4.28 -7.41 6.21 12000 1
## beta_acs_shade[35] -1.82 4.58 -9.50 5.04 12000 1
## beta_acs_shade[36] -1.76 5.70 -11.01 7.13 12000 1
## sigma 12.72 0.29 12.28 13.20 12000 1
## sigma_acs 5.06 0.94 3.55 6.46 5823 1
## sigma_shelf 3.54 1.57 1.58 5.43 4851 1compare(m.shade.species.int.acs.shelf,m.shade.species.acs.int.shelf)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.acs.int.shelf 8055.5 63.3 0.0 1 58.42 NA
## m.shade.species.int.acs.shelf 8083.3 42.3 27.8 0 56.39 17.27coeftab(m.shade.species.acs.int.shelf)## m.shade.species.acs.int.shelf
## alpha 41.5
## alpha_acs[1] -2.65
## alpha_acs[2] -4.97
## alpha_acs[3] -4
## alpha_acs[4] 3.02
## alpha_acs[5] 8.86
## alpha_acs[6] 0.18
## alpha_acs[7] 6.43
## alpha_acs[8] 0.88
## alpha_acs[9] 2.83
## alpha_acs[10] 0.39
## alpha_acs[11] 3.07
## alpha_acs[12] -5.6
## alpha_acs[13] -0.32
## alpha_acs[14] 1.88
## alpha_acs[15] -2.3
## alpha_acs[16] 2.69
## alpha_acs[17] 4.89
## alpha_acs[18] -8.27
## alpha_acs[19] 6.31
## alpha_acs[20] -3.45
## alpha_acs[21] 0.24
## alpha_acs[22] 1.56
## alpha_acs[23] -1.91
## alpha_acs[24] -2.62
## alpha_acs[25] -2.02
## alpha_acs[26] -1.83
## alpha_acs[27] 1.13
## alpha_acs[28] -0.9
## alpha_acs[29] -0.36
## alpha_acs[30] 4.25
## alpha_acs[31] 6.29
## alpha_acs[32] -0.07
## alpha_acs[33] 0.04
## alpha_acs[34] -7.22
## alpha_acs[35] -4.39
## alpha_acs[36] -2.68
## alpha_shelf[1] -3.33
## alpha_shelf[2] 4
## alpha_shelf[3] 0.13
## alpha_shelf[4] -2.4
## alpha_shelf[5] -1.34
## alpha_shelf[6] 2.54
## beta_sp[1] -3.61
## beta_sp[2] 2.37
## beta_sp[3] 3.01
## beta_sp[4] -11.41
## beta_sp[5] 6.84
## beta_shade 22.08
## beta_acs_shade[1] -13.4
## beta_acs_shade[2] -10.71
## beta_acs_shade[3] -3.19
## beta_acs_shade[4] -5.14
## beta_acs_shade[5] 3.22
## beta_acs_shade[6] 0.45
## beta_acs_shade[7] 6.54
## beta_acs_shade[8] 11.25
## beta_acs_shade[9] -1.5
## beta_acs_shade[10] 10.56
## beta_acs_shade[11] -0.61
## beta_acs_shade[12] -12.25
## beta_acs_shade[13] -2.57
## beta_acs_shade[14] -10.53
## beta_acs_shade[15] -11.97
## beta_acs_shade[16] -1.49
## beta_acs_shade[17] 3.94
## beta_acs_shade[18] -2.93
## beta_acs_shade[19] 17.64
## beta_acs_shade[20] -5.75
## beta_acs_shade[21] 0.53
## beta_acs_shade[22] 10.26
## beta_acs_shade[23] 4.65
## beta_acs_shade[24] -7.81
## beta_acs_shade[25] -10.29
## beta_acs_shade[26] 8.61
## beta_acs_shade[27] 6.76
## beta_acs_shade[28] -3.36
## beta_acs_shade[29] 8.4
## beta_acs_shade[30] 8.65
## beta_acs_shade[31] 5.16
## beta_acs_shade[32] 8.15
## beta_acs_shade[33] -0.57
## beta_acs_shade[34] -0.7
## beta_acs_shade[35] -1.82
## beta_acs_shade[36] -1.76
## sigma 12.72
## sigma_acs 5.06
## sigma_shelf 3.54
## nobs 1008plot(coeftab(m.shade.species.acs.int.shelf),pars=c(1:42))
plot(coeftab(m.shade.species.acs.int.shelf),pars=c(43:86))
We don’t need both a species response and an accession response (interaction) term. If we are only keeping one of the terms then the accession interaction model is preferred.
From this I conclude that variation is better considered at the accession level than at the species level.
Smoking deaths among doctors
In 1961 Doll and Hill sent out a questionnaire to all men on the British Medical Register inquiring about their smoking habits. Almost 70% of such men replied. Death certificates were obtained for medical practitioners and causes of death were assigned on the basis of these certificates. The breslow data set contains the person-years of observations and deaths from coronary artery disease accumulated during the first ten years of the study.
Analyse this data set to determine the posterior probability that smoking increases death by coronary artery disease, that age increases death by coronary artery disease, and that there is an interaction between age and smoking.
You can load the data set and learn about the columns using the commands below
data("breslow",package = "boot")
help("breslow",package ="boot")
breslow## age smoke n y ns
## 1 40 0 18790 2 0
## 2 50 0 10673 12 0
## 3 60 0 5710 28 0
## 4 70 0 2585 28 0
## 5 80 0 1462 31 0
## 6 40 1 52407 32 52407
## 7 50 1 43248 104 43248
## 8 60 1 28612 206 28612
## 9 70 1 12663 186 12663
## 10 80 1 5317 102 5317You can think of “person-years” as the number of observations
breslow$n <- as.integer(breslow$n)
m.smoke <- map2stan(alist(
y~dbinom(n,p),
logit(p) <- alpha + beta_smoke * smoke,
alpha ~ dnorm(0,5),
beta_smoke ~ dnorm(0,5)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 3.6e-05 seconds (Sampling)
## 4e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.smoke)
pairs(m.smoke)
precis(m.smoke)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha -5.96 0.1 -6.12 -5.81 1059 1
## beta_smoke 0.55 0.1 0.39 0.71 1071 1logistic(-5.97)## [1] 0.002547734logistic(-5.97 + 0.55)## [1] 0.004407633breslow$age_id <- coerce_index(breslow$age)
m.age <- map2stan(alist(
y~dbinom(n,p),
logit(p) <- alpha_age[age_id],
alpha_age[age_id] ~ dnorm(0,5)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 4e-05 seconds (Sampling)
## 4.4e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.age)
precis(m.age,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_age[1] -7.65 0.16 -7.92 -7.40 3234 1
## alpha_age[2] -6.14 0.09 -6.28 -6.00 3800 1
## alpha_age[3] -4.98 0.06 -5.08 -4.88 3655 1
## alpha_age[4] -4.25 0.07 -4.36 -4.15 3571 1
## alpha_age[5] -3.91 0.08 -4.03 -3.77 2974 1logistic(precis(m.age,depth=2)@output$Mean)## [1] 0.0004747746 0.0021465847 0.0068064289 0.0140105972 0.0196180444pairs(m.age)
compare(m.age,m.smoke)## WAIC pWAIC dWAIC weight SE dSE
## m.age 8614.5 4.5 0.0 1 268.92 NA
## m.smoke 9495.8 1.9 881.3 0 296.55 62.31m.smoke.age <- map2stan(
alist(y~dbinom(n,p),
logit(p) <- alpha_age[age_id] + beta_smoke*smoke,
alpha_age[age_id] ~ dnorm(0,5),
beta_smoke ~ dnorm(0,5)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 4e-05 seconds (Sampling)
## 4.5e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.smoke.age)
pairs(m.smoke.age)
precis(m.smoke.age,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_age[1] -7.92 0.19 -8.19 -7.60 1732 1
## alpha_age[2] -6.43 0.13 -6.64 -6.24 1221 1
## alpha_age[3] -5.28 0.11 -5.45 -5.10 1046 1
## alpha_age[4] -4.55 0.11 -4.74 -4.38 1058 1
## alpha_age[5] -4.20 0.12 -4.38 -3.99 1163 1
## beta_smoke 0.35 0.11 0.18 0.52 867 1compare(m.age,m.smoke,m.smoke.age)## WAIC pWAIC dWAIC weight SE dSE
## m.smoke.age 8604.7 5.6 0.0 0.99 268.47 NA
## m.age 8614.5 4.5 9.8 0.01 268.92 6.60
## m.smoke 9495.8 1.9 891.1 0.00 296.55 61.06m.smoke.age.int <- map2stan(
alist(y~dbinom(n,p),
logit(p) <- alpha_age[age_id] +
beta_smoke*smoke +
beta_smoke_age[age_id]*smoke,
alpha_age[age_id] ~ dnorm(0,5),
beta_smoke ~ dnorm(0,5),
beta_smoke_age[age_id] ~ dnorm(0,5)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 6e-06 seconds (Warm-up)
## 4.2e-05 seconds (Sampling)
## 4.8e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.smoke.age.int)
precis(m.smoke.age.int,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_age[1] -9.19 0.70 -10.24 -8.04 2301 1.00
## alpha_age[2] -6.82 0.29 -7.29 -6.39 2988 1.00
## alpha_age[3] -5.32 0.19 -5.65 -5.05 3377 1.00
## alpha_age[4] -4.53 0.19 -4.83 -4.21 2918 1.00
## alpha_age[5] -3.85 0.18 -4.14 -3.56 3210 1.00
## beta_smoke 0.53 2.04 -2.79 3.68 1004 1.00
## beta_smoke_age[1] 1.25 2.13 -1.99 4.77 1049 1.01
## beta_smoke_age[2] 0.26 2.06 -3.05 3.52 1018 1.00
## beta_smoke_age[3] -0.13 2.05 -3.34 3.14 1005 1.00
## beta_smoke_age[4] -0.21 2.05 -3.57 2.89 1000 1.00
## beta_smoke_age[5] -0.62 2.05 -4.01 2.50 1008 1.00compare(m.smoke,m.age,m.smoke.age,m.smoke.age.int)## WAIC pWAIC dWAIC weight SE dSE
## m.smoke.age.int 8601.6 10.0 0.0 0.83 268.38 NA
## m.smoke.age 8604.7 5.6 3.1 0.17 268.47 7.01
## m.age 8614.5 4.5 12.9 0.00 268.92 9.49
## m.smoke 9495.8 1.9 894.2 0.00 296.55 61.21
pred <- ensemble(m.smoke.age,m.smoke.age.int)
mu <- apply(pred$link,2,mean)
PI <- apply(pred$link,2,PI)
plot.frame <- data.frame(age=breslow$age,smoke=as.factor(breslow$smoke),mean=mu,low=PI[1,],high=PI[2,])
plot.frame## age smoke mean low high
## 1 40 0 0.0001677758 0.0000357960 0.0003888649
## 2 50 0 0.0012077086 0.0006502804 0.0017695528
## 3 60 0 0.0049718544 0.0037737892 0.0063698872
## 4 70 0 0.0108106227 0.0077777626 0.0144388881
## 5 80 0 0.0200786706 0.0136071769 0.0268564114
## 6 40 1 0.0005961793 0.0004272274 0.0007811791
## 7 50 1 0.0023859608 0.0020253972 0.0027634900
## 8 60 1 0.0072170599 0.0064127637 0.0080622621
## 9 70 1 0.0146490377 0.0129864984 0.0164390941
## 10 80 1 0.0195327976 0.0165407160 0.0225883393pl <- ggplot(plot.frame,aes(x=age,y=mean,ymax=high,ymin=low,fill=smoke))
pl <- pl + geom_bar(position = "dodge",stat = "identity")
pl + geom_errorbar(position = position_dodge(width=0.9),width=0.5)
breslow$age_num <- as.numeric(as.character(breslow$age))
m.smoke.age.int.numeric <- map2stan(
alist(y~dbinom(n,p),
logit(p) <- alpha +
beta_age * age_num +
beta_smoke*smoke +
beta_smoke_age*smoke*age_num,
alpha ~ dnorm(0,10),
c(beta_smoke,beta_age) ~ dnorm(0,5),
beta_smoke_age ~ dnorm(0,1)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 4.2e-05 seconds (Sampling)
## 4.7e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.smoke.age.int.numeric)
precis(m.smoke.age.int.numeric)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha -12.06 0.54 -12.95 -11.22 407 1.02
## beta_smoke 2.05 0.57 1.11 2.94 395 1.02
## beta_age 0.11 0.01 0.09 0.12 412 1.02
## beta_smoke_age -0.03 0.01 -0.04 -0.01 404 1.02compare(m.smoke.age.int, m.smoke.age.int.numeric)## WAIC pWAIC dWAIC weight SE dSE
## m.smoke.age.int 8601.6 10.0 0.0 1 268.38 NA
## m.smoke.age.int.numeric 8647.8 3.5 46.2 0 269.84 15.02
Cane Sugar
This data comes from an experiment to measure disease resistance in different varieties of sugar cane.
You can get the data and learn about it with:
data("cane",package="boot")
help("cane",package="boot")
head(cane)## n r x var block
## 1 87 76 19 1 A
## 2 119 8 14 2 A
## 3 94 74 9 3 A
## 4 95 11 12 4 A
## 5 134 0 12 5 A
## 6 92 0 3 6 Asummary(cane)## n r x var block
## Min. : 29.00 Min. : 0.00 Min. : 1.00 1 : 4 A:45
## 1st Qu.: 85.75 1st Qu.: 3.00 1st Qu.: 7.00 10 : 4 B:45
## Median :112.50 Median : 11.50 Median :11.50 11 : 4 C:45
## Mean :118.14 Mean : 20.26 Mean :11.94 12 : 4 D:45
## 3rd Qu.:144.50 3rd Qu.: 25.25 3rd Qu.:15.00 13 : 4
## Max. :243.00 Max. :131.00 Max. :36.00 14 : 4
## (Other):156Is there evidence of differences in disease resistance in the different varieties? Does including an adaptive prior for Block improve the model fit?
mean(cane$r/cane$n)## [1] 0.1741397logit(mean(cane$r/cane$n))## [1] -1.556568cane$var_id <- coerce_index(cane$var)
cane$block_id <- coerce_index(cane$block)
cane$n <- as.integer(cane$n)
m.cane.alpha <- map2stan(
alist(
r ~ dbinom(n,p),
logit(p) <- alpha,
alpha <- dnorm(-1.55,1)
),
data=cane,
chains = 4,
cores = 2
)## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used##
## SAMPLING FOR MODEL 'r ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000141 seconds (Sampling)
## 0.000145 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.plot(m.cane.alpha)
precis(m.cane.alpha)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha -1.58 0.02 -1.6 -1.54 3562 1logistic(-1.58)## [1] 0.1707955
m.var <- map2stan(
alist(
r ~ dbinom(n,p),
logit(p) <- alpha + alpha_var[var_id],
alpha <- dnorm(-1.55,1),
alpha_var[var_id] ~ dnorm(0,5)
),
data=cane,
chains = 4,
cores = 2,
iter=10000,
warmup = 1000
)## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used##
## SAMPLING FOR MODEL 'r ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000165 seconds (Sampling)
## 0.000169 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 3600 / 36000 ]
[ 7200 / 36000 ]
[ 10800 / 36000 ]
[ 14400 / 36000 ]
[ 18000 / 36000 ]
[ 21600 / 36000 ]
[ 25200 / 36000 ]
[ 28800 / 36000 ]
[ 32400 / 36000 ]
[ 36000 / 36000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.var,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
precis(m.var,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha -1.90 0.58 -2.88 -1.01 1072 1
## alpha_var[1] 3.47 0.59 2.56 4.47 1142 1
## alpha_var[2] 0.14 0.59 -0.79 1.12 1126 1
## alpha_var[3] -2.87 0.82 -4.14 -1.54 2117 1
## alpha_var[4] 0.99 0.59 0.06 1.97 1128 1
## alpha_var[5] 1.17 0.59 0.24 2.13 1106 1
## alpha_var[6] 1.07 0.58 0.17 2.05 1088 1
## alpha_var[7] -0.71 0.61 -1.66 0.30 1194 1
## alpha_var[8] 0.43 0.58 -0.50 1.39 1102 1
## alpha_var[9] 0.12 0.59 -0.81 1.08 1110 1
## alpha_var[10] -3.26 0.83 -4.55 -1.92 2057 1
## alpha_var[11] 0.68 0.59 -0.24 1.65 1116 1
## alpha_var[12] -0.12 0.59 -1.01 0.89 1129 1
## alpha_var[13] 1.18 0.58 0.27 2.16 1087 1
## alpha_var[14] -0.05 0.59 -0.99 0.91 1119 1
## alpha_var[15] 1.72 0.58 0.85 2.72 1090 1
## alpha_var[16] -2.58 0.77 -3.80 -1.33 1835 1
## alpha_var[17] -0.25 0.59 -1.19 0.72 1135 1
## alpha_var[18] -1.16 0.63 -2.15 -0.13 1268 1
## alpha_var[19] 0.11 0.59 -0.82 1.11 1142 1
## alpha_var[20] -1.93 0.64 -2.94 -0.88 1331 1
## alpha_var[21] -0.66 0.60 -1.60 0.33 1152 1
## alpha_var[22] -1.04 0.65 -2.09 -0.01 1333 1
## alpha_var[23] 3.22 0.59 2.30 4.20 1132 1
## alpha_var[24] -0.03 0.60 -0.97 0.96 1144 1
## alpha_var[25] -4.75 1.20 -6.59 -2.87 3332 1
## alpha_var[26] -1.12 0.60 -2.07 -0.14 1151 1
## alpha_var[27] 0.38 0.59 -0.55 1.37 1139 1
## alpha_var[28] -0.05 0.59 -0.96 0.94 1119 1
## alpha_var[29] -2.05 0.68 -3.14 -0.97 1433 1
## alpha_var[30] -1.57 0.64 -2.61 -0.56 1297 1
## alpha_var[31] 0.96 0.59 0.00 1.90 1129 1
## alpha_var[32] 0.00 0.59 -0.96 0.94 1111 1
## alpha_var[33] -0.37 0.60 -1.33 0.60 1146 1
## alpha_var[34] 0.43 0.59 -0.47 1.43 1116 1
## alpha_var[35] 1.49 0.58 0.55 2.42 1091 1
## alpha_var[36] -1.24 0.65 -2.34 -0.24 1368 1
## alpha_var[37] -0.30 0.59 -1.22 0.69 1132 1
## alpha_var[38] -1.26 0.63 -2.29 -0.27 1250 1
## alpha_var[39] -2.19 0.66 -3.21 -1.10 1371 1
## alpha_var[40] 1.29 0.58 0.40 2.28 1105 1
## alpha_var[41] -1.85 0.64 -2.87 -0.83 1273 1
## alpha_var[42] 0.81 0.58 -0.11 1.78 1103 1
## alpha_var[43] -0.25 0.60 -1.19 0.73 1161 1
## alpha_var[44] 0.94 0.59 0.00 1.90 1113 1
## alpha_var[45] 1.41 0.59 0.49 2.39 1112 1compare(m.cane.alpha,m.var)## WAIC pWAIC dWAIC weight SE dSE
## m.var 15759.3 45.8 0.0 1 170.57 NA
## m.cane.alpha 19489.1 1.0 3729.9 0 173.20 115.5
m.var.block <- map2stan(
alist(
r ~ dbinom(n,p),
logit(p) <- alpha + alpha_var[var_id] + alpha_block[block_id],
alpha <- dnorm(-1.55,1),
alpha_var[var_id] ~ dnorm(0,5),
alpha_block[block_id] ~ dnorm(0,sigma_block),
sigma_block ~ dcauchy(0,1)
),
data=cane,
chains = 4,
cores = 2,
iter=10000,
warmup=4000
)## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used##
## SAMPLING FOR MODEL 'r ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.000176 seconds (Sampling)
## 0.000179 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 2400 / 24000 ]
[ 4800 / 24000 ]
[ 7200 / 24000 ]
[ 9600 / 24000 ]
[ 12000 / 24000 ]
[ 14400 / 24000 ]
[ 16800 / 24000 ]
[ 19200 / 24000 ]
[ 21600 / 24000 ]
[ 24000 / 24000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")